function [handles] = ProcessImageGUI(varargin) %ProcessImageGUI(handles)

% [handles] = ProcessImageGUI(handles)
% Top-level image analysis function for our experiment.
% Input:
% handles - A structure that must contain at least the following values:
% handles.do2DFit - boolean
% handles.dpath - string, full path to file
% handles.autoROI - boolean
% handles.do1DFit - boolean
% handles.ROI - [1 x 4] matrix of integers
% handles.useBackground - boolean
% handles.scaleBackground - boolean
% If it contains handles.figure1, it is assumed the function is
% called from the GUI and will act accordingly (try to update GUI
% display).
% Returns:
% handles - Same as input structure but with following values added:
% handles.numAtoms
% handles.pX
% handles.pY
% Note 'handles' can just be a path to the HDF file, in which case the
% function assumes default values for the handles structure.

% TODO:
%   Save to file:
%       Integrals of x- and y- data (but over which ROI?)
%   Choose fitting function
%   if 'handles.do1DFit' is 0, then don't display flat lines

if ischar(varargin{1})
    handles.dpath = varargin{1};
    handles.do2DFit = 0;
    handles.autoROI = 0;
    handles.do1DFit = 1;
    handles.ROI = [1 1392 1 1040];
    handles.useBackground = 0;
    handles.scaleBackground = 0;
    handles.saveOD = 0;
    doDisplayImages = 1;
elseif isstruct(varargin{1})
    handles = varargin{1};
else
    error('Input must be a string or structure')
end

% To check if function has been called from GUI or not,
% if it has, update the GUI appropriately, if not, just
% analyse the data and pass the results out in the structure.
% TODO: save analysis results in HDF file?
calledFromGUI = isfield(handles,'figure1');

% Inputs from hdf5 file
[rawimages,img_time,type_img,expp,pxsize] = ReadImage(handles.dpath);
% TODO: CHANGE LABVIEW SO THAT CCDbin IS SAVED IN HDF5 FILE
%       MODIFY ReadImage()ACCORDINGLY

% Libraries used in this function
cons = Constants();

[ODraw,ovflw] = TreatImage(rawimages,type_img,handles);

% Filter image to find ROI
ODfilt = FilterImage(ODraw); % Low-pass filter image

yC = 1:expp.short_px_nb;   %HORIZONTAL DIMENSION OF EXPERIMENT
xC = 1:expp.long_px_nb;

[xROI,yROI] = BuildROI(ODfilt,expp,handles);

yInt = sum(ODraw(yROI,xROI),2);  %sum along x-direction
xInt = sum(ODraw(yROI,xROI),1);

if handles.do1DFit % Real fit of ODraw using initial guesses
    % INTEGRATED ODraw in diff. dimens. using autoROI as size (for real fit)
    pYtmp(1:3) = [max(yInt)-min(yInt) mean(yROI) 150];
    pXtmp(1:3) = [max(xInt)-min(xInt) mean(xROI) 150];
    
    options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
    pY = fminsearch(@(p)FitGauss(p,yROI,yInt),[pYtmp(1:3) min(yInt)],options);
    pX = fminsearch(@(p)FitGauss(p,xROI,xInt'),[pXtmp(1:3) min(xInt)],options);
else
    pY = zeros(1,4);
    pX = pY;
end

if handles.do2DFit
    [ODrawtreated,bindim]=TreatODfor2Dfit(ODraw,xROI,yROI);
    
    axes(handles.axes_ODraw_4_analysis)
    cax = PlotSurf(ODrawtreated);
    
    [IC2dfit,IniCond2d] = IniConditions2DFit(ODrawtreated); % Using ODrawtreated for 2D fitting IC
    
    [fitcoordmatrix,xdatabin,ydatabin,Z2dbin] = For_nlinfit2D(ODrawtreated); % Creating matrices for 2D fitting
    
    % Fitting - data is treated with equal weight
    opts=statset('TolFun',1e-6);
    [P2dbin,r2dbin,J2dbin]=nlinfit(fitcoordmatrix,Z2dbin(:),@gauss2derror,IniCond2d,opts);
    P2d = [P2dbin(1),P2dbin(2),P2dbin(3)*bindim(2),P2dbin(4)*bindim(1),P2dbin(5)*bindim(2),P2dbin(6)*bindim(1),P2dbin(7)/bindim(2),P2dbin(8)/bindim(1)];
    
    % Generate OD plot based on fit parameters
    [plotcoordmatrix,xdata,ydata,Z2data] = For_nlinfit2D(ODraw(yROI,xROI)); % y-dir. is HORIZ. exper. dir.; x-dir. is VERT
    
    ODcalc = ODcalculate2D(P2d,plotcoordmatrix,xdata);  %OD calculated from fit paramters
    ODcalcbin = ODcalculate2D(P2dbin,fitcoordmatrix,xdatabin);  %"Binned" OD calculated from fit paramters
    ODcalcIC = ODcalculate2D(IniCond2d,fitcoordmatrix,xdatabin);  %OD calculated from initial conditions
    ODresiduals = ODraw(yROI,xROI)-ODcalc;
    
    axes(handles.axes_ODcalc_2Dfit)
    PlotSurf(ODcalc,cax);  %Plot 2d fit using fit param
    
    axes(handles.axes_ODresid_rawmcalc)
    PlotSurf(ODresiduals,cax);  %Plot 2d fit residuals
else
    P2d = [0,0,0,0,0,0,0,0];
end

handles = AnalysisToHandles(pX,pY,P2d,pxsize,type_img,img_time,rawimages,ODraw,xROI,yROI,expp,cons,handles);

SaveImgAnalysis(type_img,handles,rawimages,yInt,xInt,ODraw)

if calledFromGUI % If GUI called this function then update GUI display
    axes(handles.axes_processed)
    imagesc(ODraw)
    ROIlines(xROI,yROI)
    axis equal
    xlims = xlim;
    ylims = ylim;
    set(gca,'FontSize',10)
    
    axes(handles.axes_y_proj)
    plot(yInt,yROI,gauss(pY,yC),yC);
    set(gca,'YDir','reverse');
    ylim(ylims);
    set(gca,'FontSize',10)
    
    axes(handles.axes_x_proj)
    plot(xROI,xInt,xC,gauss(pX,xC));
    xlim(xlims);
    set(gca,'FontSize',10)
    
    DisplayRawImages(type_img,rawimages,handles)
    
    DispTextOnGUI(type_img,pX,pY,handles,xROI,yROI,cons)
    
elseif doDisplayImages
    figure(1)
    subplot(3,3,[1 2 4 5])
    imagesc(ODraw)
    ROIlines(xROI,yROI)
    axis equal
    xlims = xlim;
    ylims = ylim;
    set(gca,'FontSize',10)
    
    subplot(3,3,[3 6])
    plot(yInt,yROI,gauss(pY,yC),yC);
    set(gca,'YDir','reverse');
    ylim(ylims);
    set(gca,'FontSize',10)
    
    subplot(3,3,[7 8])
    plot(xROI,xInt,xC,gauss(pX,xC));
    xlim(xlims);
    set(gca,'FontSize',10)
end % end calledFromGUI